syms q qdot qdotdot l mass g real

i = [1 0 0]';
j = [0 1 0]';
k = [0 0 1]';

% g = magnitude of gravity
G = -g*j;

F = mass*G;

% Location of the center of mass relative to the newtonian frame.
center = [-l*q 1 0]';

% Location of the point mass expressed in newtonian frame
rp = center + l*[cos(q) sin(q) 0]';

% Velocity of the point mass
vp = diff(rp,q)*qdot;

% Acceleration of the point mass
ap1 = diff(vp,q)*qdot;
ap2 = diff(vp,qdot)*qdotdot;
ap =  ap1 + ap2;

% Partial Velocity
vp_r = diff(vp,qdot);
 
Fr = dot(F,vp_r)-dot(mass*ap1,vp_r);
lhs = diff(dot(mass*ap2,vp_r),qdotdot); 
udot = Fr/lhs;

% Expression for KE and PE of the system
ke = 0.5*mass*(vp'*vp);
pe = mass*g*rp(2);